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; SUMMARY 

i 

A process for the compatibility check of measured aircraft responses is 
presented. This process includes estimation of bias errors in the data and 
comparison of reconstructed responses with those measured directly. The model 
relating aircraft states and outputs is based on the six-degree-of-freedom 
kinematic equations and on output equations specifying the measured variables. 
Three model forms are presented. One of them includes the effect of atmo- 
spheric turbulence. 

i ■ 

The estimation technique used is an extended Kalman filter with one-stage 
optimal smoothing. The algorithm also includes a time-varying process-noise 
covariance matrix, a fixed-point smoother, and an analysis of residuals. 

The resulting technique is applied to simulated data to check the accuracy 
of the computing algorithm and the estimated bias errors. The compatibility 
check is also applied to measured flight data obtained from longitudinal and 
lateral maneuvers of a general aviation aircraft. Improved accuracy of measured 
data is obtained when the data are corrected for estimated bias errors. 


INTRODUCTION 

Estimation of aerodynamic characteristics of aircraft from flight data is 
receiving increased attention, partially because of the availability of sophisti- 
cated algorithms and computers for this purpose. In practice, it is found that 
biases or unknown scale factors are usually present in recorded data. Erroneous 
results for the aerodynamic characteristics are obtained and valuable research 
and computer time may be wasted before the researcher finds that the various 
data channels are not compatible. For example, data incompatibility would exist 
if the measured incidence angles did not agree with those reconstructed from the 
accelerometer and rate gyro measurements. A more efficient method for flight- 
test data analysis would be to determine first whether the data channels are all 
compatible, make whatever adjustments are required, and then use the refined 
"compatible data" for parameter estimation. 

One of the first rigorous attempts of checking data compatibility was made 
in reference 1 , in which the main objective was flight-path reconstruction from 
test data. Regression analysis and redundancy in measured data from three- 
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degree-of-freedom longitudinal motion were used. This work was followed by ref- 
erence 2, in which the aircraft states, constant biases in measured data, and 
initial conditions were estimated. The results were obtained from maximum like- 
lihood estimation and an extended Kalman filter combined with fixed-interval 
smoothing. Later, using the latter technique, an extensive study of the same 
problem was made in reference 3. 

A practical technique for obtaining the weighted least squares estimates of 
aircraft states, measurement biases, and initial conditions from measured data 
was introduced in reference 4. The mathematical model used was formed by the 
six-degree-of-freedom kinematic equations with measured accelerations and angu- 
lar rates as input variables. The same model was used in reference 5, but the 
estimation algorithm was comprised of an extended Kalman filter-smoother and a 
fixed-point smoother. This technique was developed in reference 6 as "an 
advanced method" for identification of aircraft parameters and states. 

The present report describes a computing procedure which 

( 1 ) Checks the consistency of the measured responses of an aircraft after 

removing the estimated bias errors in the measured data 

(2) Reproduces some of the state variables and their first derivatives 

which are not measured directly 

(3) Estimates the means, variances, autocorrelation functions, and power 

spectral densities of residuals 

The model for the system to be analyzed is based on the aircraft six-degree-of- 
freedom kinematic equations and output equations specifying the measured varia- 
bles. Three different models are considered. One of them includes the effect 
of atmospheric turbulence. This generalization has not been previously reported 
in the literature. 

For estimation of aircraft states and bias errors, an extended Kalman fil- 
ter from reference 6 was adopted. The algorithm was modified by introducing a 
time-varying process-noise covariance matrix and an analysis of residuals. The 
first modification, which extends the method of reference 6, improves the accu- 
racy of the estimates in those cases where the process noise is nonstationary. 
The analysis of residuals should improve the filter performance by providing 
an indication for the selection of process- and measurement-noise covariance 
matrices. At the same time, the residuals give reasonable information about 
measurement-noise characteristics because the mathematical model used for 
the aircraft state estimation is, in general, very well defined. 

The technique presented here should form part of a comprehensive process 
for the analysis of flight-test data and should precede any attempts at aircraft 
parameter estimation. 
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SYMBOLS 


a x> a y> a z 

by 

C 1 > c 2 
D 

E{ } 

, f 

f( ) 
g 

i H 

i 

1 h 
h( ) 

I 

K 

| L 
m 
N 
n 

n-| ,n2, . . . 
P 

P»Q»r 

Q 
R 

R v (t) 

Sv(f) 


longitudinal, lateral, and vertical components of acceleration, m/sec 
constant bias error in variable y 

constants in differential equation for Gauss-Markoff process 

matrix of dummy variables in fixed-point smoother algorithm 

expected value 

frequency, sec"^ 

function which represents model 

acceleration due to gravity, m/sec 2 

3h( )/3x* evaluated along reference trajectory (see eq. (B2)) 
altitude, m 

nonlinear output vector used to represent measurement system 

identity matrix 

Kalman filter gain 

scale of turbulence, m 

mean value of residuals 

number of data points 

measurement-noise vector 

,n-|5 elements of measurement-noise vector 
covariance matrix of augmented state variables 
roll, pitch, and yaw velocities, rad/sec or deg/sec 
process-noise covariance matrix 
jth main diagonal element of Q matrix 
measurement-noise covariance matrix 
autocorrelation function of residuals v 
power spectral density of residuals v 


S(w) 

s 

t 


u,v,w 


u g> v g’ w g 


V 


X 

x« 

Xb»yb> z b 

x a »Ya 




y 

z 


a v 

&v 

<5(t) ,§ij 

n 

0 

e 

Xy 

v y 

K 

a 

a 2 

x 


power spectral density of turbulence components 
standard error estimate of residuals 
time, sec 

longitudinal, lateral, and vertical airspeed components, m/sec 

longitudinal, lateral, and vertical turbulence velocity components, 
m/sec 

true airspeed, m/sec 
state vector 
augmented state vector 

linear position coordinates of aircraft, m 

position coordinates of a v wind vane with respect to aircraft cente: 
of gravity, m 

position coordinates of 8 V wind vane with respect to aircraft cente 
of gravity, m 

output vector 

measurement vector 

angle of attack measured by wind vane, rad or deg 
sideslip angle measured by wind vane, rad or deg 
Kronecker delta for continuous and discrete data, respectively 
input vector 

vector of unknown parameters 

pitch angle, rad or deg 

scale factor error of variable y 

residual of variable y 

process-noise vector 

standard error 

variance 

time lag, sec 
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transition matrix for linearized equations 


tp roll angle, rad or deg 

j yaw angle, rad or deg 

w angular frequency, rad/sec 

Subscripts : 

E measured quantity 

g gust 

i index of data points 

R uncorrected for bias errors 

0 initial value 

Superscript: 

(j) index of iterations 

Matrix exponents : 

T transpose matrix 

-1 inverse matrix 

Mathematical Notation: 

j • over symbols denotes derivative with respect to time 

over symbol denotes estimated value 
A incremental value 

A symbol followed by (i|i) indicates the conditional mean or variance. 
For example, x(i-l|i) is the estimate of Xj__-| based upon all measured data 
up to time t^. 

MATHEMATICAL MODEL RELATING AIRCRAFT STATES AND OUTPUTS 

The mathematical model used for the data compatibility check is described, 
in general, by three sets of kinematic equations with the state variables con- 
sisting of three linear velocities u, v, and w; three position angles 9, 0, 

and v|»; and three linear positions X5, y^, and Z5. The input variables in 

the equations are the linear accelerations a x , a y , and a z and angular rates 
p, q, and r. The kinematic equations (state equations) are formulated as 



u = -qw + rv + a x - g sin 0 ( 1 ) 
v = -ru + pw + ay + g cos 0 sin ip (2) 
w = qu - pv + a z + g cos 6 cos <p (3) 
<p = p + q sin 9 tan 0 + r cos q> tan 0 (4) 
0 = q cos tp - r sin tp (5) 
tp = q sin tp/cos 0 + r cos tp/cos 0 (6) 


x b = u cos 0 cos ip + v(sin tp 
+ w(cos 9 sin 0 cos tp + 


sin 0 cos tp - cos tp sin tp) 
sin tp sin tp) 


Yb = u cos tp sin tp + 
+ w(cos tp sin 0 


v(sin tp sin 0 sin tp + cos tp cos tp) 
sin tp - sin tp cos tp) 


(7) 


( 8 ) 


zjj = -u sin 0 + v cos 0 sin tp + w cos 0 cos tp (9) 

It is considered that the following variables are measured: 

(1) The inputs to the system a x , ay, a z , p, q, and r 

(2) The airspeed V, two incidence angles 0 V and oi v , three position 

angles tp, 0, and tp, and altitude h = -z^ (these variables repre- 
sent the output of the system) 

The measured variables z are corrupted by systematic and random errors. It 
is, therefore, assumed that each can be expressed as 

z = (1 + Xy)y + by + n (10) 

where y is the true value of the output, Xy is the unknown scale factor, by 
is the constant bias error, and n is the measurements noise. It is further 
assumed that the scale factor Xy is nonzero only for variables V, 8 V , and 
a v which are directly connected with airflow measurement. 

In flight-test experiments, the horizontal and lateral positions of an air- 
craft with respect to a given point on the ground are usually neither measured 
nor requested. For that reason, the system of state equations can be simpli- 
fied by deleting equations (7) and (8). Then, replacing the input variables in 
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the remaining state equations by their measured values results in a new set of 
state equations as follows: 


u 


0 

r R,E - b r 

“Ur,E - b q) 

o’ 


u 

4 

V 


-(r R)E - b r ) 

0 

PR,E “ b p 

0 


V 

w 


9 R,E - b q 

-(PR,E “ b p> 

0 

0 


w 

• 

h 


sin 0 

-cos 0 sin 9 

-cos 0 COS <P 

0 


h 




-g sin 0 



a x,R,E 


~ b ax 


0 

w 

-V 

6 


n 4 


n i 


g 

cos 0 sin 

9 


a y,R,E 


_b ay 


-w 

0 

u 

0 


n 5 


n 2 

+ 

g 

cos 0 cos 

9 

+ 

a z,R,E 


-b az 

+ 

V 

-u 

0 

0 


n 6 

+ 

n 3 



0 



0 


0 


0 

0 

0 

0 


0 


0 


(11a) 


9 


1 sin cp tan 0 cos 9 tan 0 


PR,E " b p + n 4 

0 

= 

0 cos tp -sin 9 


9 R,E - b q + n 5 



0 sin tp/cos 0 cos 9 /cos 0 


r R,E - b r + n 6 


! The output equations take the form 


Vp s (1 + + by 


0 V ,R = (1 + *g) tan -1 
a v,R = (1 + tan -1 



+ bg 

+ b a 


9r = 9 + b tp 

0r s 6 + bg 

^R = ^ + 

h R = h + b h 


( 11 b) 


( 12 ) 
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where by are the constant bias errors in the input and output variables. Equa- 
tions (11) and (12) are similar to those developed and used in reference 5. 

In formulating the output equations for 8 V R and a v R , the effect of ran- 
dom errors in p, q, and r was neglected. The output equations also do not 
include the effect of atmospheric disturbances. Therefore, it is further 
assumed that all flight data are obtained from measurement in still air. How- 
ever, a gust disturbance model which can be incorporated into the analysis is 
presented in a subsequent section. 

The measured values in the state equations introduce process noise into the 
system and make the system model stochastic. The process noise in equations (11) 
is nonstationary and cross-correlated. Two assumptions are made concerning the 
process noise: 

(a) All cross-correlations between process-noise components are neglected. 

(b) The process noise has zero mean value and its covariance matrix varies 
with time. The expressions for the elements of the covariance matrix are given 
in appendix A. 

The general form of the state equation for the given system can be written 
as 


x = f(x,n,©,t) + £(t) 


(13) 


and the discrete form of the measurement equation as 


Zj_ = h i (x i ,n i ,0) + ni (i = 0,1,. . . ,N) (14) 

where x is the state vector, n is the input vector, z is the measurement 
vector, 0^ = [by,Xy] is the vector of unknown parameters, £(t) is the 
process-noise vector with E|fj(t)} = 0 and E|fj(t) £^(t)| = Q(t) 6(t-x), n^ 
is the measurement-noise vector with Ejn-jJ = 0, and E|n^nj^| = R<Sjj, and N 
is the number of data points. The matrices Q(t) and R are the covariance 
matrices for the process and measurement noise, respectively, the first one 
being time dependent and the other constant. 


Model With Measured Angular Accelerations 

In some cases, the angular accelerations p, q, and r can also be 
obtained from measurement. When substituted into the system model equations, 
these accelerations form the input variables, and the angular rates become out- 
puts of the system. The state equations are changed as follows: 
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• 

u 


0 

r 

-q 

o' 


u 



-g sin 0 

* 

V 


-r 

0 

P 

0 


V 

+ 

g 

cos 0 sin cp 

• 

w 


q 

-P 

0 

0 


w 


g 

cos 0 cos cp 

h 


sin 0 

-cos 0 sin <p 

-cos 0 COS ip 

0 


h 



0 


a x, R,E - b ax + n 1 

a y , R , E “ b ay + n 2 

a z, R,E “ b az + n 3 


(15a) 


• 

P 


PR,E - b p + n lH 


q 

= 

qR,E “ b q + n 15 


» 

r 


^R,E - bp + n-,6 




1 sin <p tan 0 

cos q> tan 0 

0 

s 

0 cos (p 

-sin ip 



0 sin cp/cos 0 

cos cp/cos 0 


(15b) 


(15c) 


where b£, b^, and bf are the constant bias errors in the angular accelera- 

tions. Output equations (12) are completed by three equations: 

PR = P + b p 


q R = q + bq 


r R = r + t>r 

In equations (15), the process noise enters explicitly and can be consid- 
ered stationary. Additional measurements of angular accelerations provide possi- 
bilities for estimating the measurement-noise characteristics of angular veloci- 
ties and also the process-noise characteristics, as explained in reference 6. 
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Model With Gust Disturbances 


In general, gusts can be taken into account by modeling the spectra of atmo- 
spheric turbulence. The model is based on the assumption that the turbulence 
field is steady and homogeneous. By using Dryden's model (see, e.g., ref. 7, 
pp. 315-318), the power spectral densities of the horizontal, lateral, and verti- 
cal turbulence velocity components are represented by 


^ug ) 


Syg (W ) 


S W g(w) 



( 16 ) 


(17) 


( 18 ) 


where L u , L v , and L w are the scales of turbulence and a U g, a vg> an ^ 0 
are the variances of the turbulence velocity components. It is also assumed 
that the turbulence field is isotropic; therefore, 


L u = 


L w = L 


2 2 2 
a ug = a vg = a wg = a g 2 

For further simplification of the turbulence modeling, equations (17) and (18) 
can be approximated by 


Syg (W ) 


S w g(w) 


°s 2 


3L 

2tt 


1 


1 + 


L_ 

2 


to \ £ 

V 


(19) 


Turbulence with the power spectral density given by equation (16) or (19) 
can be described as a Gauss-Markoff process of first order. The corresponding 
differential equation is given as 

Xg = -c-jXg + C2ng (20) 


where x g is a turbulence velocity component u g , v g , or w g and n g repre- 
sents a normally distributed white noise with zero mean and variance The 
constants c-) and 03 have the form for horizontal gust, 
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°2 = 



and the lateral and vertical gust, 



The only measurements which the gust affects directly are airspeed, side- 
slip angle, and angle of attack. The perturbations due to gust on these varia- 
bles may be expressed as 




'g u g 






OLr 


w, 


g 


( 21 ) 


With the gust disturbances taken into account, state equations (11) or (15) 
can be extended by three equations for u g , v g , and w g based on equation (20). 
The first three output equations can be modified to 


Vr s (1 + Xy) (v/ u 2 + v^ + 




^R 


(1 + Xg) tan- 1 1 v * (r B>E ~ b r> x P ~ ( PR,E - 

(1 + X a ) tan" 1 j " ~ ^R,E “ bq)xq + (p R)E - V/q 



+ bg 



+ b a 


For the solution of the extended system of state equations, the values of L 
and O g 2 must be given, A good approximation for the scale of turbulence is 
L = 300. A rough estimate of <J g 2 can be obtained from the random component 
of the recorded wind-vane reading. 


IDENTIFICATION TECHNIQUE 


By using the extended Kalman filter, identification of the system given by 
equations ( 13 ) and ( 1 4) can be achieved as the state estimation of a nonlinear 
system with the augmented state vector 


x* = 



( 22 ) 
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Because of augmenting the state vector with the unknown parameter vector, state 
equation (13) is changed to 


f(x,q ,0 ,t) 

+ 

S(t)" 

0 


. o . 


(23) 


The algorithm is obtained by linearizing the system and measurement 
equations around the best estimate of states at each data point. In order to 
correct for system and measurement nonlinearities, the computing algorithm 
presented in reference 5 includes "local iteration" and one-stage optimal 
smoothing. At time t^_i , the algorithm starts with x*( j) (i-1 1 i-1 ) and 
p( j) ( i_1 | i_1 ) , where the exponent j indicates the jth iteration. Then the 
algorithm predicts and updates states and state variances to time t^ by apply- 
ing one iteration of the extended filter. Finally, on the basis of this new 
estimate, it smooths back to t^_i . The smoothing closes the loop to provide 
an improved reference for the new prediction at t^. The iteration terminates 
when there is no significant difference between consecutive iterations or after 
a specific number of iterations. Analysis in reference 6 showed that this proce- 
dure can significantly reduce the bias which is inherent in the extended Kalman 
filter. The computing scheme for one stage with two iterations is presented in 
figure 1 . 

At each data point the extended Kalman filter provides the information for 
the fixed-point smoother. The basic feature of the algorithm is that no storage 
for the filtered states and state covariance matrices is required. The algo- 
rithm works in conjunction with the extended Kalman filter. After each step 
from t^_i to t^, the smoother computes the smoothed estimates x*(o|i) and 
P ( 0 1 i ) , thus ending with x*(0|N) and P ( 0 1 N ) when all data are processed. 

The algorithm for the locally iterated extended Kalman filter combined with 
the fixed-point smoother is summarized in appendix B. The schematic block dia- 
gram for the identification technique indicating initial data required and 
results obtained is given in figure 2. 


ANALYSIS OF RESIDUALS 


Residuals are defined as the differences between measured variables and 
their predicted values: 

v i = z i ~ y(i| i_1 ) (24) 

For the correct values of the process-noise and measurement-noise covariance 
matrices and the correct model of the system, the residuals should approach a 
random Gaussian, white sequence with zero mean. The variances should be consis- 
tent with calculated values from filter equation (B5); that is, 


E 


I (j) (j) T 
W Vi 




pU)(i|i-l) 


H 


(j) T 


+ R 


(25) 
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Because of its presence in equation (B4), the process-noise covariance 
matrix influences the predicted state covariance matrix. Therefore, equa- 
tion (25) can be used for the consistency check between the matrices and 

R and the time histories of residuals, provided that the model defined by 
equations (23) and (14) is correct. The plotted residuals should be within the 
boundaries given by the square root of the main diagonal elements in the matrix 
defined by equation (25). If the boundaries are either too large or too small, 
the matrices and R which are not known a priori and must therefore be 

estimated in advance, should be adjusted accordingly. 

The mathematical model of the system is based on the general form of the 
aircraft kinematic equations. The uncertainty in the model used is in the form 
of the measurement equation (10), where only constant bias errors are consid- 
ered. Despite this uncertainty, the assumption of negligible modeling errors 
can be substantiated. Then, with good estimates for and R, the residuals 

reflect the characteristics of the measurement-noise components, mainly the fre- 
quency contents in their power spectral densities. 

For these reasons mentioned, the extended Kalman filter-smoother has been 
complemented by the analysis of residuals. From their time histories, the mean 
values, variances, autocorrelation functions, and power spectral densities of 
the residuals are estimated. For the estimates of the autocorrelation functions 
and power spectral densities, the expressions from reference 8 (pp. 290-295) 
were used. 


EXAMPLES USING SIMULATED DATA 

! 

The data consistency check procedure developed in the present study was 
r applied to simulated data to check the accuracy of the computing algorithm and 
the estimated bias errors. Simulated data used represent the six-degree-of- 
freedom maneuvers of an aircraft. A summary of the measured variables and 
standard errors of the simulated measurement noise is given in table I. In 
table II, the elements of the augmented state vector and the initial estimates 
of their variances are presented. The number of data points is N = 400 and 
the time interval is At = 0.05 sec. 

First the simulated data were used to evaluate the effect of varying the 
number of terms in the transition matrix (eq. ( B1 ) ) and the number of iterations 
in the single-stage smoother. As a result, six terms of the Taylor series expan- 
sion were used for the approximation of the state transition matrix, and the num- 
ber of iterations in the single-stage smoother was set equal to two. 

The effect of using a time-dependent or constant process-noise covariance 
matrix was also investigated. The initial estimates 

©0 T = [ b y» X y] 

and x*(o| 0) were set equal to their true values, and two cases were computed. 
Results are given in table III. In case 1 the Q matrix was kept constant; 
that is, Qi = Qo for* all i* Ih case 2 the Q matrix varied with time as 
expressed by equations (A3). Comparison of results shows improved accuracy of 
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the parameter estimates in case 2, although their standard errors (lower bounds) 
presented in parentheses are larger than those in case 1. It was therefore 
decided to include equations (A3) in the algorithm. Similar procedures in ref- 
erences 2, 3, and 4 use a constant approximation for the Q matrix. 

In table IV the effect of unknown initial values for the parameter esti- 
mates and their variances and the number of unknown parameters is demonstrated. 
For case 3, the data were processed five times. The initial estimates of the 
unknown parameters were first set equal to zero. Then, in the following four 
passes, they were always updated with the final estimates from the previous 
pass. For several parameters, improved accuracy of the estimates after five 
passes is apparent. This example also shows the small effect of the initial 
values P ( 0 | 0 ) on the final estimates P(n|n). In case 4 the number of unknown 
parameters was reduced from 14 to 6. The new set of unknown parameters was con- 
nected with variables defining the lateral motion of an aircraft. The initial 
parameter estimates were again set equal to zero. The decreased number of 
unknowns resulted in better accuracy of the estimates, although it was not possi 
ble to estimate the parameter Xg more accurately. This was caused by the low 
sensitivity of the output variables to the parameter Xg as indicated by its 
standard error. 


ANALYSIS OF MEASURED FLIGHT DATA 

Flight-test data were obtained on a small single-engine general aviation 
aircraft. The first run analyzed represents longitudinal motion of the aircraft 
excited by elevator deflection. In the second example, measurements were avail- 
able from eight runs using rudder and aileron deflections as input. For these 
runs, flight conditions were the same; only the forms of both inputs were differ 
ent. The use of large amplitudes of rudder and aileron deflection resulted in 
coupling between the lateral and the longitudinal motion of the aircraft. 


Three-Degree-of-Freedom Longitudinal Motion 

The model for the extended Kalman filter-smoother included the following 
state, input, and output variables: 

x T = [u,w,6] 

= [ a x,R,E> a z,R,E><)R,E] 
y T = [ v R> a v,R* S R] 

The vector of unknown parameters was assumed in the form 
= j^bg X , b a 2 » bq , by , b(j , bQ , Xy , X^J 

The time histories of measured input and output variables are plotted in fig- 
ure 3- For the measured data, the sampling interval was At = 0.1 sec and the 
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number of data points was N * 198. The estimated parameters, standard errors, 
and sensitivities from the three oomputer passes are given in table V. In the 
first pass, the initial values of the parameters were set equal to zero and the 
initial variances were set equal to those given in table VI. Then, in the fol- 
lowing passes, the initial values were always equal to the final fixed-point 
smoother estimates from the preceding pass. Results from the first and third 
passes show some changes in parameters and slightly improved accuracy of the 
final estimates. The accuracy is expressed as the Cramer-Rao lower bound on 
the standard error of the estimates. These values are given in parentheses. 

Some of the estimated parameters have high standard error which can be due 
to their small significance in the model. In order to investigate this possi- 
bility, the sensitivities of the output variables with respect to the unknown 
parameters were computed. The expression for the sensitivity was taken from 

reference 9 as where Pjj is the main diagonal element of the infor- 

mation matrix P _1 (N|N). The computed sensitivities are also presented in 

table V. 

In the fourth pass the parameters b az , b q , and be with high uncertainty 
and small sensitivity were fixed at their last estimated values. The resulting 
estimates of remaining parameters were more accurate and had more uniform sensi- 
tivities than those from the complete model. 

For all passes the initial values of the process-noise covariance matrix 

Qq were the same. The initial values of the measurement-noise covariance 

matrix were changed by comparing their values with the estimates P(o|N) for 

the state variables and with the variances of the residuals. 

The time histories of the parameters from the fourth pass are shown in fig- 
ure 4. With the exception of b ax , they do not indicate any substantial changes 
toward the end of the time interval. The corresponding variances were almost 
constant within the whole time interval. 

In table VI the initial values of the state variables and variances from 
two passes are compared with the smoothed values. There is only slightly better 
agreement between these values after two passes. 

The measured and predicted output time histories are plotted in figure 5 , 
and the residuals in figure 6. The 2a bounds in the residual plots were esti- 
mated from equation (25). The agreement between measured and reconstructed 
responses indicates good data compatibility. Because almost all residuals for 
each variable are within 20 bounds, reasonable initial values for the matrices 
R and Qq can be expected. The outputs from the deterministic model based on 
equations (1) to (5) were also computed and compared with measurements in fig- 
ure 7. The resulting errors are given in figure 8. The errors in the determin- 
istic model are quite large in all variables, especially in a v . The main rea- 
son might be due to neutrally stable modes in the deterministic model. 

The analysis of residuals is summarized in figures 9 and 10, where the mean 
values, standard errors, autocorrelation functions, and power spectral densities 
of residuals are presented. From figures 9 and 10, it is apparent that only the 
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residuals in the pitch angle are close to the random white sequence. The resid- 
uals in the remaining two variables are distorted by deterministic components 
of an unknown origin. 


Six-Degree-of-Freedom Motion 

The compatibility check was also applied to measurements from eight runs in 
which large-amplitude rudder and aileron deflection resulted in coupling between 
the lateral and longitudinal motion of the aircraft. The sampling interval for 
the measured data was At = 0.05 sec; the number of measured data points was 
N = 200 for the first four runs and N = 1400 for the remaining runs. The 
state, input, and output vectors were defined as 

x T = [u,v,w,(p,6] 

r 1 T = £ a x , R , E > a y , R , E > a z , R , E > PE > <3E > r E J 

y T = [Vr.Bv.R^v.R/PR* 0 ]*] 

In the first pass the vector of unknown parameters had 12 elements. After the 
uncertainty and sensitivity of all the estimated parameters were examined, the 
number of unknowns was reduced by two and the vector © was formed as 

©^ = [^>ax > b ay » ^az > ^r * > tyx > fyp , b© , ^y , 

The estimated parameters from eight runs are presented in table VII. 

The standard errors of the same parameters differed little between runs. There- 
fore, only the average values of these standard errors from all eight runs are 
included. For some parameters, the estimates are quite consistent. The largest 
inconsistency was observed in the bias error for the sideslip angle. 

The comparison of the measured and predicted outputs for run 5 is given 
in figure 11. The comparison of measured and reconstructed variables is good. 
Similar agreements were observed in the remaining runs. The residuals from 
run 5 are plotted in figure 12, and the autocorrelation functions in figure 13. 
As for three-degree-of- freedom longitudinal motion, the residuals in some vari- 
ables indicate the presence of deterministic components. 


CONCLUDING REMARKS 

A process for the compatibility check of measured aircraft responses has 
been presented. This process includes the estimation of bias errors in the data 
and comparison of some reconstructed responses with those measured directly. 

The model relating aircraft states and outputs is based on the six-degree-of- 
freedom kinematic equations and on output equations specifying the measured 
variables. Three model forms are presented, one of which includes the effect 
of atmospheric turbulence. 
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The estimation technique used was an extended Kalman filter with one-stage 
optimal smoothing. Application of this smoothing can significantly reduce the 
bias which is inherent in the extended Kalman filter. In order to further 
improve the accuracy of the estimates, the algorithm includes a time-varying 
process-noise covariance matrix, a fixed-point smoother, and an analysis of 
residuals. The smoother gives smoothed initial estimates of states, parameters, 
and their covariance matrix. The analysis of residuals can help to verify the 
preselected values of the process- and measurement-noise covariance matrices by 
comparing the residual variances obtained from the filter with those from the 
time histories of residuals. In addition, it provides estimates of characteris- 
tics of the measurement noise in terms of the autocorrelation functions or power 
spectral densities. 

The resulting technique was first applied to simulated data to check the 
accuracy of the computing algorithm and the estimated bias errors. It was shown 
that the replacement of a constant process-noise covariance matrix with the time- 
variable one can improve the estimates where the process noise is nonstationary. 

The accuracy of the estimates was also improved by iteratively processing 
the data several times. The reprocessing used the final estimates from the 
previous pass as the new initial estimates for parameters. 

The compatibility check was also applied to measured flight data. In the 
first example the data front longitudinal motion were analyzed. The covariance 
matrix of the states and parameters and the inverse of this matrix, were used 
for the assessment of the relative accuracy of the estimates and the sensitivity 
of the model outputs with respect to unknown parameters. Then, in the repeated 
pass, parameters with very small accuracy and sensitivity were dropped from the 
model. This simplification of the model resulted in improved accuracy of the 
remaining parameters. 

The second example included a compatibility check of the lateral responses 
of an aircraft in eight repeated runs. The repeatability of estimates obtained 
was good even when model included 5 states and 10 parameters. 

The estimated autocorrelation functions and power spectral densities of 
residuals in both examples show that in some cases the estimated measurement- 
noise characteristics may be different from those assumed for white noise. The 
reason for this discrepancy might be any remaining unknown deterministic compo- 
nents in the measured data. 

w#; ' 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, VA 23665 
Jime 14, 1977 
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APPENDIX A 


PROCESS-NOISE COVARIANCE MATRIX 
The process-noise vector in state equation (23) has the form 

*11 
^2 
0 


where 


•1 = 


0 

-w 

v 

0 


w 

0 


u 

0 

0 


r 

n 4 


V 


n 5 

+ 

n 2 

i 

n 6 


n 3 

1 

j 

0 _ 


_0 


( A 1 ) 


^2 


COS (p 


-sin tp 


1 

n 4 


n 5 

I 

_ n 6_ 


(A2) 


By using equations (A1) and (A2) and the assumption of no correlation 
between the measurement-noise components, the main diagonal elements of the 
process-noise covariance matrix can be expressed as 


<111 


( 1 2 + w 2 dg 2 


v 2 a g 2 






(A3) 


<122 = <^2^ + w 2 a 4 2 + u 2 dg 2 

<133 = a 3 2 + + u 2 dg 2 

944 = 0 

qgg = a 4 2 + (sin cp tan 0) 2 ag 2 + ( cos tan 0)2ag 2 
qgg = cos 2 ip CJ5 2 + sin 2 q> dg 2 
qjy = (sin ip/cos 0)2a5 2 + (cos cp/cos 0)2a g 2 

<188 = Q99 = . . . = 923 23 = 0 

J 

where dj 2 is the variance of the measurement noise nj for j = 1,2,. . .,6. 
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Equations (A3) can be simplified by replacing the time variables 
u,v,. . . ,(cos (p/cos 0) by their average values which depend upon the flight 
record being analyzed. 
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APPENDIX B 


EXTENDED KALMAN FILTER-SMOOTHER AND FIXED-POINT SMOOTHER 

In the algorithm from reference 5 for an extended Kalman filter-smoother, 
the following notation is introduced: 


$(j>(i,i-1) 


I + §£ At + At 2 + 

9x 2\9x*/ 


J !x*( J) (i-1 | i) 


H - ( j ) = x*( J)(i| i)] 

1 " 57 * 

ViU) = - hj_ [&*( J) (i | i)J - Hj/ jjx*^ J) (i| i-1 ) - x*(j)(i|i)J 


(B1 ) 
(B2) 


At = tj_ - t^_i 

Before the computing starts, it is necessary to set i = 1, j = 1, and 
x*( 1 ) (i-1 | i) = x#(i-l|i-1). Then the filter-smoother proceeds in the following 
steps: 

(1) State prediction 

x*( j) (i| i-1 ) = f ' f(x*,ri,t) dt + j)(i,i-1 ) [x*(i-1 1 i-1 ) - x*( j) ( i-1 1 i)] 


Ai.-i (B3) 

If j = 1, set x*^j)(i|i) = x*(J)(i|i-1) 

(2) State covariance prediction 

p(J)(i|i-1) = *<J)(i,i-1) P ( i— 1 | i-1 ) j) T (i , i-1 ) + Qj. (B4) 

(3) Gain estimation 

Kj/J) = p(J)(i|i-1) Hi^ J) T |^H i ( j) p(J)(i|i-1) Hi^J) 1 + r] (B5) 

(4) State updating 

x*( j +1 ) ( i | i) = x*< J) (i| i-1 ) + K i U)v i (j) (B6) 

(5) State covariance updating 

P<J)(i|i) = (i - K i (j)H i (J))p(j)(i|i-l) (B7) 
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(6) State smoothing 

x«( j+1)(i_i |i) = x«(i-1 |i-1 ) + P(i-1 ji-1) $(j) T (i,i-i) 


[l^U) p(j)(i|i-1) HiU) 1 + r] 


-1 


Vi (J) 


(B8) 


When the first iteration is completed, j is incremented and the filter- 
smoother starts again with equation (B3). If x*(J)(i|i) = x*(J +1 )(i|i) or 
j = J ma y which is the specified number of iterations, then x*(i|i) is set 
equal to x*(J +1 )(i|i) and P(i|i) equal to p(j)(i|i), and the filter- 
smoother continues with the incremented value for i. 


At the last iteration of the filter-smoother, the values of %, v if 
<&(i,i-1), and are used in the fixed-point smoother. In order to initialize 

this part of the algorithm, the matrix Dq with dummy variables is set equal to 
P(0|0). Then 

Di = Di_i $ T (i,i-1) (I - KiHi) 1 (B9) 

x*(0| i) = x*(0| i-1 ) + D i H i T R~ 1 V i (BIO) 


P(0|i) = P(0|i-1) - DiHjTR- 1 <D(i,i-1) HiDi.-, 


(BID 
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TABLE I.- STANDARD ERRORS OF SIMULATED MEASUREMENT NOISE 


Variable 

Standard error o 
of measurement 
noise of variable 

a x , m/sec 2 . . . 

0.02 

a y , m/sec 2 . . . 

0.02 

a z , m/sec 2 . . . 

0.10 

p, rad/sec ... 

0.0008 

q, rad/sec . . . 

0.0008 

r, rad/sec . . . 

0.0008 

V, m/sec .... 

0.20 

! 

B v , rad .... 

0.0002 

a v , rad .... 

0.0008 

q> , rad . . . . . 

0.0008 

0, rad ..... 

0.0008 







TABLE II.- ELEMENTS OF AUGMENTED STATE VECTOR AND INITIAL ESTIMATES 


OF AUGMENTED STATE VARIABLE VARIANCES FOR SIMULATED DATA 


Augmented state variable 

i 

Initial estimate of 

variance of variable 
(a) 

u, m/sec 

1 .0 

v, m/sec 

1.0 

w, m/sec 

1 .0 

ip, rad 

0.00002 

9, rad 

0.00002 

b ax , m/sec 2 

0.10 

b a y, m/sec 2 

0.10 

b az , m/sec 2 

0.10 

bp, rad/sec 

0.00005 

bq, rad/sec 

0.00005 

b r , rad/sec 

0.00005 

by, m/sec 

0.10 

bp , rad 

0.0001 

b a , rad 

0.0001 

bqj , rad 

0.0001 

bp , rad 

0.0001 

Xy 

0.00001 



0.00001 

x a 

0.00001 


a Main diagonal elements of P ( 0 | 0 ) . 



TABLE III.- EFFECT OF PROCESS-NOISE COVARIANCE MATRIX FORM 


ON PARAMETER ESTIMATES FOR SIMULATED DATA 


Parameter 

True value 

Estimate of parameter 
(a) 

Case 1 

Q = Constant 

Case 2 
Q = Q(t) 

} m/ S6C^ • • • 

0.20 

0.665 (0.0032) 

0.25 (0.045) 

b ay , m/sec^ . . . 

0.20 

0.215 (0.0038) 

0.30 (0.15) 

b az , m/sec^ . . . 

1.0 

0.505 (0.0045) 

0.96 (0.15) 

bp, rad/sec . . . 

0.004 

0.00407 (0.000018) 

0.0040 (0.0010) 

bq, rad/sec . . . 

0.004 

0.00327 (0.000022) 

0.0033 (0.00089) 

b r , rad/sec . . . 

0.004 

0.00428 (0.000014) 

0.0042 (0.00090) 

by, m/sec .... 

2.0 

1.75 (0.057) 

1.89 (0.29) 

bg, rad 

0.002 

-0.01620 (0.000084) 

0.0084 (0.0030) 

b a , rad 

0.01 

0.00628 (0.00010) 

0.010 (0.0032) 

fyp > rad 

0.01 

0.0143 (0.00013) 

0.0090 (0.0032) 

be , rad 

0.01 

0.00849 (0.000085) 

0.010 (0.0026) 

x v 

0.10 

0.1027 (0.00049) 

0.10 (0.0027) 

x 3 

0.10 

0.0306 (0.00032) 

0.097 (0.0031) 

x a 

0.10 

0.1107 (0.00022) 

0.10 (0.0028) 


a Numbers in parentheses are Cramer-Rao lower bounds on standard 
errors. 
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TABLE IV.- EFFECT OF INITIAL VALUES OF UNKNOWN PARAMETERS AND 


REDUCED NUMBER OF UNKNOWNS ON PARAMETER ESTIMATES 
FOR SIMULATED DATA 


Parameter 

True 

value 

Estimate of parameter 
(a) 

Case 3 

Case 4 

1st computer 
pass 

5th computer 
pass 

j m/ see • • • 

0.20 

-0.36 (0.047) 

0.081 (0.046) 


bay » m/sec 2 . . . 

0.20 

0.78 (0.16) 

0.67 (0.15) 

0.23 (0.18) 

b az , m/sec 2 . . . 

1 .0 

1.42 (0.16) 

1.14 (0.15) 


bp, rad/sec . . . 

0.004 

0.0059 (0.0010) 

0.0046 (0.0010) 

0.0038 (0.00099) 

bq, rad/sec . . . 

0.004 

0.0030 (0.00087) 

0.0032 (0.00088) 


b r , rad/sec . . . 

0.004 

0.0044 (0.00089) 

0.0047 (0.00090) 

0.0035 (0.00096) 

by, m/sec .... 

2.0 

1.74 (0.29) 

3.10 (0.29) 


bg, rad 

0.002 

0.0024 (0.0027) 

0.0049 (0.0027) 

0.0052 (0.0029) 

b/M . rad 

0.01 

0.022 (0 0029) 

0011 (0 0031 ) 


b^>, rad 

0.01 

0.0060 (0.003D 

0.0078 (0.0032) 

0.0094 (0.0032) 

bA . rad 

0.01 

0.0065 (0.0026) 

0.0097 (0 0026) 


Xy 

0. 1 

0.034 (0.0026) 

0.084 (0.0027) 




0.1 

0.0003 (0.0031) 

0.0011 ( 0 . 0031 ) 

0.0002 (0.0031) 


0. 1 

0.027 (0.0027) 

0.073 (0.0027) 







a Numbers in parentheses are Cramer-Rao lower bounds on standard errors. 
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TABLE V.- ESTIMATES AND SENSITIVITIES OF PARAMETERS 
FROM FOUR COMPUTER PASSES FOR FLIGHT DATA 




1st computer pass 

2nd computer pass 

4th computer pass 


Parameter 

Estimate 

of 

parameter 

(a) 

Estimate 

of 

parameter 

(a) 

Sensitivity 

of 

parameter 

Estimate 

of 

parameter 

(a) 

Sensitivity 

of 

parameter 

b ax> 

m/sec 2 . . . 

-0.063 

(0.031) 

-0.063 

(0.028) 

5.0 

-0.065 

(0.019) 

12.3 

h „ 

m/sec 2 . . . 

0.0012 

(0.081) 

0.0016 

(0.079) 

0.0004 



u az> 




rad/sec . . . 

-0.00007 

(0.00088) 

-0.00007 

(0.00086) 

0.006 



Dq> 



by, 

m/sec .... 

-0.060 

(0.087) 

-0.076 

(0.080) 

1.5 

-0.134 

(0.076) 

8.2 

ba* 

rad 

0.0023 

(0.0024) 

0.0024 

(0.0023) 

1.7 

0.0022 

(0.0021) 

1.8 

b 0» 

rad .... 

0.0003 

(0.0022) 

0.0004 

(0.0021) 

0.07 






x v 


0.0031 

(0.0022) 

0.0019 

(0.0021) 

1.7 

0.0047 

(0.0018) 

18.2 

X a 


0.0019 

(0.003D 

1 

0.0021 

(0.0031) 

0.43 

0.0041 

(0.0031) 

1.7 


a Numbers in parentheses are Cramer-Rao lower bounds on standard errors. 
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TABLE VI.- INITIAL AND SMOOTHED VALUES OF STATE VARIABLES AND 
VARIANCES FROM THREE COMPUTER PASSES FOR FLIGHT DATA 


State variable 

Pass 

no. 

Estimate of variable 

Variance of variable 

Initial value, 
x(0 | 0) 

Smoothed value, 

x(o|n) 

Initial value 
(a) 

Smoothed value 
(b) 

u, m/sec . . . 

1 

46.174 

46.09 

0.01 

0.0076 

w, m/sec . . . 

1 

2.368 

2.37 

0.01 

0.0083 

0 , rad .... 

1 

0.03941 

0.0394 

0.000001 

0.00000095 

u, m/sec . . . 

3 

46.174 

46.16 

0.0076 

0.0064 

w, m/sec . . . 

3 

2.368 

2.37 

0.0083 

0.0080 

0 , rad .... 

1 3 

0.03941 

0.0394 

0.00000095 

1 i 

0.00000094 


a Main diagonal elements of P ( 0 | 0 ) . 
b Main diagonal elements of P ( 0 | N ) . 
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Figure 1.- Computing scheme for two iterations of extended 
Kalman filter-smoother. 
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Figure 2.- Block diagram 
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Figure 5.- Comparison of measured and predicted time histories of output 
variables. Three-degree-of-freedom longitudinal motion. 
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Figure 7.- Measured and computed time histories of output variables from 
deterministic model. Three-degree-of- freedom longitudinal motion. 
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Time, t, sec 

Figure 8.- Time histories of errors in output variables from deterministic 
model. Three-degree-of-freedom longitudinal motion. 
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Figure 9.- Mean values, standard errors, and autocorrelation functions of 
residuals. Three-degree-of- freedom longitudinal motion. 
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Figure 10.- Power spectral densities of residuals. Three-degree-of-freedom 

longitudinal motion. 
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Figure 11.- Comparison of measured and predicted time histories of output 
variables for six-degree-of-freedom motion. Flight-data run 5. 
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Figure 11.- Concluded. 
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